Phyllotaxis, Pushed Pattern-Forming Fronts and Optimal Packing 
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We demonstrate that the pattern forming partial differential equation derived from the auxin 
distribution model proposed by Meyerowitz, Traas and others gives rise to all spiral phyllotaxis 
properties observed on plants. We show how the advancing pushed pattern front chooses spiral 
families enumerated by Fibonacci sequences with all attendant self similar properties and connect 
the results with the optimal packing based algorithms previously used to explain phyllotaxis. Our 
results allow us to make experimentally testable predictions. 
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INTRODUCTION 

Using a model derived from the pioneering ideas of 
Meyerowitz, Traas et al [1] (MT) on the role of PIN 1 
proteins in creating an instability of uniform auxin con- 
centrations near a plant's shoot apical meristem (SAM), 
this Letter reports on stunning new results which lend 
credence to the view that almost all of the features of 
phyllotactic configurations are the result of a pushed pat- 
tern forming front whose origin is the MT instability. The 
front leaves in its wake either whorls or Fibonacci spirals. 
The patterns we observe exhibit all known self-similar 
properties, reveal some new invariants and reproduce the 
well known van Iterson diagrams associated with discrete 
algorithms (Levitov, Douady and Couder, Atela, Gole 
and Hotton Q) which reflect optimal packing strategies 
based on ideas and observations of Hofmeister and Snow 
and Snow Q. Further, the location of the maxima of 
the auxin fields coincide very closely with the point con- 
figurations generated by the discrete algorithms, which 
suggests that pattern forming systems may be a new tool 
for addressing optimal packing challenges. In short, in- 
stability generated patterns are the mechanism by which 
plants and other organisms can pursue optimal strate- 
gies. 

The equation we use [4] for the auxin concentration 
fluctuation about its mean derives from the continuum 
approximation to a discrete cell by cell description pro- 
posed in The resulting PDE is strikingly similar to 
those found in many pattern forming systems which sug- 
gests that Fibonacci spiral patterns should be observable 
in many physical contexts. The instability of the uniform 
auxin concentration state which gives rise to the pattern 
is primarily due to reverse diffusion. In most equilib- 
rium situations, inhomogeneities in chemical concentra- 
tions are smoothed by ordinary diffusion. But in plants 
the situation is not an equilibrium one. PIN 1 proteins in 
the interiors of cells move under the influence of an auxin 
gradient to the cell walls where they orient so as to drive 
auxin in the direction of its gradient. When the effect 
is sufficiently strong, the net diffusion is negative and an 
instability with a preferred length scale occurs. As the 




FIG. 1. A pseudocolor plot of u(x.,t) on the inner sec- 
tion r < 89 of a pattern initiated at r = 233 with 
parastichy numbers (89,144). A movie may be found at 
http: //math. arizona. edu/~ pennybacker /medi a/ sunf lower/ 



shapes initiated by the instability grow, nonlinear inter- 
actions determine the preferred planforms. The resulting 
PDE for the auxin fluctuation concentration w(x, t) turns 
out to be very close to a gradient flow, and is given by 
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where x = (x, y) or (r, 0) are horizontal coordinates on 
the tunica (plant skin) surface, t is time, the most lin- 
early unstable wavelength is 27r, and /i and f3 are the 
coefficients of the linear growth and quadratic terms re- 
spectively. 



2 



METHODS AND DATA 

We study solutions of (pQ) in many geometries defined 
by surfaces of revolution r = r(z), r'{z) = on a cylin- 
der and r'{z) = 1 on a disc, but in this Letter we focus 
on what we call the sunflower situation in a disc geom- 
etry. Sunflowers are formed in two stages. In the first, 
flowers are initiated in a generative annulus surround- 
ing the plant's SAM and, as the plant grows, the gen- 
erative annulus and the region of phylla move further 
out, and their configurations evolve into spiral families. 
At a certain point, however, the center region consist- 
ing of mushy undifferentiated cells begins to undergo a 
phase transition and, from the outside in, create florets 
or seeds. In this process, the seeds are laid down annulus 
by annulus by an inwardly moving front and the partic- 
ular pattern which forms at a given radius stays at that 
radius. This means that any optimal packing property 
which the pattern manifests when it is first laid down re- 
mains visible. The diameter of the plant when the process 
is completed (over a time scale of several days) is of the 
order of millimeters. Thereafter the plant grows adiabati- 
cally until it reaches its mature diameter of 10-15cms. To 
simulate this situation, we initiate a spiral pattern with 
parastichies (m, n) in a circle of radius m + n for values 
of (m, n) ranging from (55, 89) to (89, 144) and allow the 
pattern to propagate inwards according to (pp) (See Fig. 
[T]). We then analyze the resulting field ^(x, t). 

The numerical algorithm used is nontrivial and is in- 
formed by the pioneering recent works of Furihata and 
Matsuo [5|. The key idea is to use the gradient property 
of (pQ), by rewriting ([2]) directly in discrete form, and 
choosing the expression for the discrete gradient which 
preserves the main dissipation property dE/dt < and 
ensures both stability and accuracy. For most choices of 
a discrete approximation for the continuum expression 
for 5£/5u, this essential property would be compromised. 
Details are given in [6|. 

Guided by previous analytical results at near-onset 
conditions, we understood several key points. First is the 
fact that Fibonacci spirals, and indeed whorl structures, 
are very much a consequence of the presence of a sign 
reversal asymmetry and the fact that the pattern is laid 
down by a moving front annulus by annulus in an envi- 
ronment of slowly changing metric. The asymmetry gives 
hexagonal lattices a special role in planar patterns be- 
cause modes exp(ik m -x), having k m = &o, with wavevec- 
tors 120 degrees apart reinforce each other via quadratic 
interactions. Likewise, in a circular geometry, for pat- 
terns laid down annulus by annulus, we observe triads of 
modes having wavevectors k m = (ra,£ m ), k n = (n,£ n ), 
k m+n = k m + k n with shapes exp(—i j £ m dr — im6) : a 
good approximation for r large, m/r finite. The signs 
of the radial wavenumbers alternate. These modes can 
reinforce each other via quadratic interactions at certain 
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FIG. 2. The maximum values of the amplitude of all inte- 
ger modes 8 < m < 89. Note that all Fibonacci amplitude 
maxima are equal, and their 2 nd harmonics and "irregular" 
sequence amplitudes are much smaller. 
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FIG. 3. The invariant amplitude curve for ampli- 
tudes a mj with {mj} = {1,2,3,5,8,13,21,...}. The 
horizontal scaling emphasizes the self-similar prop- 
erty a rnj (r) — a m]+1 (^r). A movie may be found at 



http : //math . arizona . edu/~pennybacker /media/ amplitude/ 



radii where they have almost equal amplitudes and en- 
ergetically preferred wavenumbers, with lengths given by 
km = tfn + m 2 /r 2 . But, as the radius increases, the 
lengths of the wavevectors begin to deviate from their 
optimal values. Moreover, the energetically preferred 
radial wavenumbers (f m (r),4(0) evolve along a com- 
putable locus. As a consequence, the first three of the 
new quadratically generated modes, k2 m , k2 n , k2 m + n 
move away from, whereas the last wavevector, k m+ 2 n , 
moves towards the critical circle. This results in the am- 
plitudes of the first three modes being slaved to, and 
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much less than, the amplitude of the last (see Fig. [2j). 
Thus, in situations where the pattern is laid down an- 
nulus by annulus, the quadratic interactions select from 
the set of all integers only that subset obeying Fibonacci 
rules. 

These considerations suggest that the quantities we 
should monitor are the amplitudes and radial wavenum- 
bers of the signal ^(x, t) = 5^ ■ u mj (r, t) exp(— zraj#)+(*), 
{raj} C Z, which we obtain by writing u rnj (r, t) as 
a mj (r, t) ex.-p(i(j) rnj ) . At any given r and t, the amplitudes 
are given by the sequence {a mj } and the radial wavenum- 
bers £ mj by d(j) mj / 'dr. We also monitor the front speed 

which is not constant but varies in a periodic fashion 
over distances separated by (/?, the golden number. Be- 
cause the chosen pattern combines several modes, one has 
to ask if they can all propagate as a synchronous front. 
There are two kinds of fronts [7[ . The first kind are pulled 
fronts whose properties are determined by conditions in 
the unstable state ahead of the front. The second cate- 
gory, to which Fibonacci spirals belong, are pushed fronts 
whose speeds and steepnesses exceed those of the pulled 
front and whose characteristics are determined by con- 
ditions behind the front. This is essential for Fibonacci 
patterns because it ensures both that the modes par- 
ticipating in the pattern structure all propagate in syn- 
chrony and introduces the important ingredient of bias 
into pattern choice. Namely the choice of pattern emerg- 
ing in the present generative annulus is influenced by the 
pattern already laid down in the previous one. This is 
important because the free energy landscape has many 
minima. The minimum which is chosen is the one for 
which the previous state of the system lies in its basin 
of attraction. The minima corresponding to Fibonacci 
spirals are metastable states, continuous in r. They are 
not the deepest minima in the landscape. Whereas the 
metastable Fibonacci pattern state is long lived, it will 
eventually tunnel its way into a lower energy minimum 
if one is available on the landscape. And eventually of 
course, the pattern will occupy a region large enough so 
that it will think it is on a plane. Very slowly, it will de- 
velop defects which will lead to patches of pure hexago- 
nal planforms mediated by defects. But coarsening takes 
long time. Moreover, in plants there may be a mecha- 
nism which roots the angular position of a primordium 
once it begins to grow. 

RESULTS 

Simulation data. In Fig. [TJ we show the simulated 
sunflower head. The transitions between pairs of paras- 
tichy numbers are smooth. For a wide range of /i, /?, the 
result is the same; we use /i = lxl0 -3 and f3 = 3 for all 
of the data that follow. The pattern makes its transitions 
from one almost hexagonal state to another via interme- 
diate two- mode structures, whose fields look rhombic. In 
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FIG. 4. The rise p and divergence angle 5 given by the local 
lattice structure at each radius. The shaded lines are the 
van Iterson tree, with selected parastichy pairs indicated. 
Inset is detail of the data for small p. A movie may be found at 
http: //math. arizona. edu/~pennybacker/media/divergence, 
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FIG. 5. The front speed is, the local energy e, and the local 
packing efficiency 77. The vertical axis has been rescaled to 
indicate relative variation from the mean value of each of these 
quantities. 

Fig. [2j we plot the maximum amplitudes (over the sun- 
flower head) of all integer modes 8 < ra < 89. The Fi- 
bonacci modes {ra^} = {8, 13, ... , 89} are dominant, al- 
though the 2 nd harmonics {2m j} and irregular Fibonacci 
modes {rrij-2 + m j}i passively slaved, have significant 
albeit small amplitudes. They measure the small differ- 
ence between a "perfect" and a pattern formed Fibonacci 
lattice. The amplitudes of all other integer modes are 
insignificant. The choice of the "winning" subsequence 
of dominant modes only depends on the initial condition. 
Different Fibonacci-like sequences of dominant wavenum- 
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bers even inherit the same self-similar properties dis- 
cussed in the following remarks. We also check in every 
case the self-similar property that a mj+1 (r(f) = a m .(r), 
and that the radial wavenumbers satisfy a similar rela- 
tion e mj+1 (rcp) = -£ mj (r), where <p = (1 + VE)/2 is the 
golden number. 

The field amplitudes. In Fig. [3j we show the new 
"amplitude invariant" which has no parallel in the dis- 
crete algorithms which only generate point configura- 
tions. We plot, as functions of r/tpi on a logarith- 
mic scale, the locus of all amplitudes {a rnj }, {rnj} = 
{3, 5, 8, . . . , 144}, as r decreases from 233 to 8. At 
r « 144, the configuration is almost hexagonal with 0144, 
agg, &55 occupying the positions £?, C, D respectively. 
As r decreases, the amplitude (2144 will very quickly de- 
crease leaving and dominant (and almost un- 
changed) for a short interval until the amplitude 
rises quickly to form the next hexagonal configuration at 
r « 89. As r continues to decrease, the pattern repeats 
with <289, &55, &34 and e&2i replacing 0144, asg, ^55 and 
(234. The graph of the locus of {a m .} in three dimensions 
(a m . +1 , a mj . , a mj ._ 1 ) is given in [6[ and shows a common 
homoclinic orbit joining the origin to itself with six legs 
on which there are two-mode dominated transitions from 
hexagons from circumferential wavenumbers raj+i, rrij, 
rrij-i to raj, raj_i, raj_2- 

T/ie lattice configurations. The overall arrangement 
of phylla on the sunflower head is not a fixed spiral lat- 
tice, but a slowly varying one. A description of this vari- 
ation may be found by assigning to each radius a "local" 
lattice. On the circle of radius r, we compute the am- 
plitudes and wavevectors of the dominant modes, which 
give the orientations of the parastichies and the shapes of 
the phylla. In fact, we can measure from the simulation 
data the set of wavenumbers {£ m } for all integers ra. The 
curves defined, as r varies, by (£ m ^n) for all consecutive 
pairs, which belong to some Fibonacci subsequence, of 
dominant circumferential wavenumbers ra, n are identi- 
cal. In Fig. 31 we choose to present one manifestation of 
this property, which allows a direct comparison with the 
van Iterson diagram, by performing a trivial change of co- 
ordinates onto the cylinder of radius r, taking the radial 
coordinate to the axial coordinate. This preserves the dif- 
ferential structure in the neighborhood of the circle. Sup- 
posing that the two dominant modes have circumferential 
wavenumbers that are coprime, the resulting lattice has 
maxima only at distinct and regularly spaced axial po- 
sitions. The axial distance p, normalized by the radius 
of the cylinder, and the angular difference S between two 
consecutive maxima, often called rise and divergence an- 
gle, uniquely identify the local lattice. This data is over- 
laid on a tree diagram, due originally to van Iterson, that 
identifies lattices which are rhombic, having at least four 
nearest neighbors of each point equidistant. Rhombic 
lattices are optimally packed, in that, for a given rise or 
divergence angle, they have the largest fraction of their 



area covered by identical circles centered at the lattice 
points. Hexagonal lattices are a special case, having six 
equidistant neighbors. It is clear that each of our local 
lattices is nearly rhombic, and moreover, that the lim- 
iting value of the divergence angle is the golden angle 
2tt/(^ 2 « 2.4. 

Self-similar quantities. In Fig. [5j we show the vari- 
ation of three important quantities that are unchanged 
by the transformation r — »> rep. The velocity v is the 
rate at which the front propagates inward. The local en- 
ergy e is the energy density of the local lattice, computed 
with the amplitudes and wavevectors measured at each 
radius. The local packing efficiency r\ is the largest frac- 
tion of the local cylinder, representing the pattern in a 
narrow annulus at radius r, that can be covered by identi- 
cal circles centered at maxima of the local lattice. When 
the pattern is nearly hexagonal, the local energy is low 
and both the front velocity and local packing efficiency 
are high. When the pattern is two-mode dominated, but 
less hexagonal, the opposite is true. Recent experimental 
work [8[ on masking developing seed patterns in sunflow- 
ers suggests that it may be possible to measure v and 
check that it varies as predicted. 

Whorl structures with wavevectors ra), (0, 2ra) 
such as the decussate ra = 2 only exist in finite r inter- 
vals. In p, we show how a decussate makes a transition 
to a spiral structure via an Eckhaus-like instability. We 
examine a range of such transitions. Finally, we remark 
that our work suggests that Fibonacci patterns are uni- 
versal, long lived but not infinitely long lived, in pattern 
forming planar systems for which the plane is tiled annu- 
lus by annulus. More generally, it will be interesting to 
explore the nature of patterns laid down by a front near 
regions in which the metric changes slowly. 

This work was supported by NSF grant DMS 0202440. 
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